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SUMMARY 


An implicit finite-difference code has been developed which solves the 
transport equations for the turbulence kinetic energy and its dissipation rate in 
generalized coordinates in three dimensions. The finite-difference equations are 
solved using the Beam-Warming algorithm. The kinetic energy-dissipation code, KEM, 
provides the closure; i.e., the turbulent viscosity for calculation of either com- 
pressible or incompressible flows. Turbulent internal flow over a backward-facing 
step has been calculated using the present code in conjunction with the Incompres- 
sible Navier-Stokes Code, INS3D. The results are in good agreement with experiments 
and two-dimensional computations of other researchers. 


INTRODUCTION 


There is a growing need for a two-equation turbulence model in generalized 
coordinates in three dimensions which can be used in conjunction with various flow 
codes that use the eddy viscosity calculated by an algebraic model of turbulence. 
Depending on the complexity of the geometry of application and that of the flow 
field, these codes should be able to choose between an algebraic and a two-equation 
turbulence model, the latter of which is described here. 

Since the transport equations for the turbulence kinetic energy, k, and its 
dissipation rate, e, are not strongly coupled with the flow-field equations (e.g., 
see refs. 1 and 2), they are considered here as a candidate for the system of turbu- 
lence equations to be solved independently of the flow-field equations. These 
equations are solved implicitly as a two-by-two system in the present study, but are 
uncoupled from the flow-field equations. This paper gives the details of such a 
model run in conjunction with the incompressible Navier-Stokes code, INS3D (ref. 3), 
for a standard test problem (internal flow over a backward-facing step in a two- 
dimensional channel). The results are compared with results from experiments 
(ref. 4) and an earlier computation (ref. 1). Discussion of detailed computations 
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(two dimensional) are given in reference 5. The model can be used with codes that 
solve either unsteady compressible or unsteady incompressible equations in general- 
ized coordinates with viscous terms in all the three coordinate directions. A minor 
modification is needed to change from incompressible to compressible formulations as 
well as from fully viscous to thin layer formulations. 


GOVERNING EQUATIONS 


The turbulence equations for k and e used here are cast in the conservation 
form. Using Einstein's index notation and the summation convention, transport 
equations for k and e can be written in Cartesian coordinates as the following 
two-by-two system, with k = (1/2)u.u. and e = (p/p)u. ,u. 

11 1 , J 1 , J 
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Th e ki netic energy production term, P, is due to the mean shear and is given by 
-pu^u.U i where U^ is the mean velocity and u^ is the fluctuating velocity. 

On application of the gradient diffusion hypothesis (Boussinesq, 1877), we can write 
the production term as 
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where 6.^ . is the Kronecker delta and the t^bulent viscosity y fc is given by 

Kolmogorov's hypothesis (ref. 6), y fc = C'pk ' Since at high Reynolds numbers the 

dissipation rate, e, can be assumed to be proportional to k^'^/a, we can write 

y, = C pK/e where C is a constant, 
t v y 


Also, y = y 


role of effective exchange or diffusion coefficients for 


and y = y + (y./o ), where 


y t /o k and y fc /a 


. ^ play the 

k and”e, respectively, and 

where a. and o can be interpreted as turbulent Prandtl numbers for the k and e 

transport processes. Hence, o k and o can be assumed to be close to unity in 

accordance wj,th expectations; the turbulent viscosity is redefined as 

y. = C f (pk /e)Re . The functions f 5 and f above take into account the low 
t y y ® y 


Reynolds number dependence of the constants C 2 and C , respectively (refs. 7 
and 8), and they are defined as follows: y 


f 2 = [1.0 - exp(-R 2 )] 2 
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The terms -2yk/d in the k-equation and (-2ye/d )f^ in the e-equation take 
into account the low Reynolds number effect near the walls and are due to refer- 
ence 9* The function f^ is given by 



where V is the friction velocity and d is the normal distance to the wall. The 
values of the various "constants" are chosen as follows: 


o k = 1.0 

C 1 = 1.44 

C 2 = 2.0 


3 



C 3 = 0.5 
C = 0.09 

V 

The coefficient of viscosity p is normalized by velocities are normal- 
ized by U m , distances by a characteristic length L, and the density p by p^. 
This results in the reference Reynolds number definition as Re^ = p^U^L/y^, normal- 
ization of k by U^, and that of e by lr/L. 

The transport equation for turbulence kinetic energy is obtained by taking the 
moment of the instantaneous form of the Navier-Stokes equations, adding and then 
averaging them (the square root of the turbulence kinetic energy represents a veloc- 
ity scale of the turbulent motion). The assumption of the gradient-diffusion model 
for k is introduced, i.e., the viscous diffusion of flux of k is proportional to 
the gradient of k. The k equation takes into account the history and transport 
effects of turbulence which the mixing length hypothesis cannot account for, since 
it is based on the assumption of local equilibrium everywhere in the flow field. 

The local equilibrium assumption is that the turbulence produced at a given point in 
the flow field is dissipated at the same rate as it is produced, thereby exerting no 
influence on the flow field at other points at a later time. Also from the mixing 
length hypothesis, the turbulent viscosity vanishes at those points in the flow 
field where the velocity gradients vanish, which is not true in general. 

Some researchers (e.g., ref. 10) have derived a transport equation for the 
length scale of the turbulent motion represented by a, and it has been used with 
the k equation in modeling a turbulent flow field (e.g., refs. 11 and 12). The 
turbulence length scale equation also accounts for the history and transport effects 
of turbulence. However, the specification of an appropriate length scale in the 
k-a system at the boundaries is a difficult task. The transport equation for the 
turbulence kinetic energy dissipation rate is derived based on physical considera- 
tions, and the number of the constants in the k-e model is one less than in the 
k-a model or any other two-equation model. The wall boundary conditions for the 
dissipation rate e are relatively easy to prescribe. Although an exact equation 
for the dissipation rate, e, can be derived from the instantaneous form of the 
Navier-Stokes equations (ref. 13), the transport equation for e used here is based 
on heuristic grounds (ref. 14) since the exact equation requires drastic model 
assumptions that yield an e equation with a highly empirical character. 

The equations for k and e given here are applicable for both the high and low 
turbulent Reynolds number, R = (pk /yejRe^. Near the wall, the high-Reynolds 
number form of the k-e system is modified as follows: (1) the viscous diffusion 
of k and e is included in the governing equations, and (2) the constants C 2 and 
C are allowed to be functions of the turbulent Reynolds number, R^, and appro- 
priate terms are added to the high-Reynolds number form of the k and e equations 
(refs. 7-9). Since the total dissipation rate is not zero at the wall and the 
isotropic dissipation there vanishes, a term for the dissipation at the wall is 
chosen which corresponds to the nonisotropic part of the energy dissipation. 
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Accordingly, the terra -2y(k/d ) is added to the right-hand side of the k equa- 
tion (ref. 9). With both the high- and low-Reynolds number forms thus built into 
the k-e system, the governing equations can be solved subject to the boundary 
conditions applied either at the wall, or away from the wall through the use of wall 
functions (ref. 15). 

Using the generalized coordinate transformation, 

T E t 


Xj * Xj(x.,t) 


where Xj = (S,n,c) and x^ = (x,y,z) the k-e system transforms to 
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The Jacobian, J is given by 


!^i 
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and the contravariant velocities are given by 


Q. = 3 , x . + U . 3 x- 

1 t*i j Xj A l 


( 2 ) 


BOUNDARY CONDITIONS 


There are two ways to incorporate the boundary conditions on k and e at the 
walls in the k-e system. In the first case, Dirichlet wall boundary conditions 
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are chosen so that k = 0 and e = 0 at the walls (the nonisotropic part of the 
dissipation rate at the wall is included in the k equation and a similar "wall" 
dissipation term is added to the e equation that permits setting e equal to zero 
at the walls (ref. 9); in the second case, the wall function approach is adopted 
which precludes having to integrate right down to the walls and thereby the 
stiffness problem associated with the first approach is avoided. In the present 
code, both options are available. Dirichlet inflow conditions and Neumann outflow 
and far-field conditions are used to complete the necessary boundary condition 
procedure. 

A wall function approach (ref. 15) is used to connect the outer region to the 
viscous sublayer. The assumption here is that the resultant tangential velocity at 
a point P, which is located near the wall outside the viscous sublayer, follows the 
logarithmic law of the wall, and that the turbulence is in local equilibrium at such 
a point, i.e., that the rate of production of the turbulence kinetic energy equals 
the rate at which it is dissipated. Then, at such a point, P, we have 



where E is a parameter which accounts for surface roughness or friction, or phe- 
nomena such as pressure gradient or mass injection through the walls, where k is 
the von Karman constant (= 0.42), Vp is the resultant ta ngen tial velocity at P, 

V is the resultant friction velocity and is equal to A /p , and where dp is the 
normal distance between P and the wall. The kinetic energy and the dissipation 
rate at P are then given by 


NUMERICAL SCHEME 


/Z 

V 
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T 
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The numerical scheme used to solve the k-e system, equation (2), is the 
implicit, noniterative approximate factorization algorithm of Beam and Warming 
(ref. 16). Before finite differencing the equations, the flux vector F y in 
equation (1) is redefined as 1 
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so that the second term is lumped explicitly with the source term, S. 
finite-difference expression for equation (2), we have 


Writing the 





Ax 

,n 


S 


(3) 


where the time-differencing is first-order (Euler) implicit, and the spatial deriva- 
tives are approximated by the central-differencing operator 6 

x i 

Linearizing the flux vector F\ locally in time, equation (3) can be written as 




(D n + AxS) - ID n - AxJ n+1 


k[('I ■ i V) • 7 Ml} 


where the Jacobian matrix of the flux vector is given by 


M. = ( A,B,C) = 3 n F n = — 
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with 



The right-hand side of equation (4) is calculated explicitly, and the left-hand side 
can be factored approximately. Therefore, we have 
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+ e^A^ (D n+1 - D n ) = RHS - e e [(V Xi A Xi ) 2 ]D n 


where 7 and A are first-order forward- and backward-difference operators, RHS is 
the right-hand side of equation (4), is the implicit second-order smoothing 

parameter, and e is the explicit fourth-order smoothing parameter. The smoothing 
terms are added to damp the high frequency oscillations in the solution which arise 
out of the central-difference approximation of the spatial derivatives. 

CODE VALIDATION RESULTS 


The k-e turbulence model code, KEM, which is developed here, is tested and 
validated by computing the flow over a backward-facing step, by comparing the pre- 
dictions with experiments (ref. 4), and by using earlier computations (ref. 1). For 
detailed simulations, the reader is referred to reference 5. 

Since the present computations were carried out with the primary purpose of 
validating the KEM code, only preliminary results showing a comparison of the pres- 
ent predictions with experimental data and with an earlier computation are given 
here. 

Figure 1 shows the plot of stream function contours over the back-step. The 
length of the separation bubble is somewhat underpredicted. A reasonable comparison 
is made of mean velocity at a location of 5.2 step heights downstream of the step as 
it is plotted against the cross flow coordinate in figure 2. Figure 3 shows the 
turbulence kinetic energy variation along the crossflow coordinate at 7.7 step 
heights downstream of the step. Finally, figure 4 shows the pressure distribution 
along the streamwise direction on the step side wall. Again, the comparison is 
reasonable. 


CONCLUSIONS 


The KEM code offers the capability of predicting the compressible or incompres- 
sible three-dimensional turbulent flows around arbitrary geometries in external and 
internal flow situations in conjunction with any flow solver. The code is robust 
and easy to use. The use of the code necessitates allocation of additional computer 
memory for the two scalar variables, the turbulence kinetic energy and the dissipa- 
tion rate. It predicts complex turbulent flows with large separated and secondary 
flow regions satisfactorily. For flows with strong pressure gradients and curvature 
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effects, the turbulence model employed here may need appropriate corrections to 
yield more accurate results. 


ADDENDUM 


A more realistic alternative boundary condition formulation has been used 
yielding results in better agreement with experimental data. This formulation is 
based on the near wall behavior of turbulence 1 ^ and yields the following boundary 
condition on e, 

£ wall = 4vk P /d P " e p 


This is similar to the deduced expression for ne ar the wall, although this 

expression has never been used as a wall boundary condition on e to the author's 
knowledge. Both the expressions due to Refs. 8 and 9, and Ref. 17 were used in the 
present study. The latter tends to be well behaved, especially in regions of small 
kinetic energy. This boundary condition provides a self-correcting mechanism in 
the k equation through the dissipation term appearing in the source terms for 
the k-equation. However, realistic initial profiles for k and e should be pro- 
vided, otherwise these boundary conditions destabilize the k-e system. It is best 
to start with a zero gradient wall boundary condition on e, and switch to this 
formulation at a later stage in the computations. It is interesting to note that if 
the zero gradient condition were true, then the boundary condition due to Ref. 9 
will be approximately the same as that due to Ref. 17 with wall gradient approxi- 
mated by the one-sided first order difference formula. 
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Figure 2.- Mean velocity profile comparison at x/h = 5.2. 
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Figure 3.- Comparison of turbulence kinetic energy at x/h = 7.7. 
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